home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 3 / Info_Mac_1994-01.iso / Science / Gnuplot 3.5 / demo / stat.inc < prev    next >
Text File  |  1993-11-03  |  6KB  |  259 lines
  1. #
  2. # $Id: stat.inc 3.38.2.6 1992/11/14 02:25:21 woo Exp $
  3. #
  4. # Library of Statistical Functions version 3.0
  5. #
  6. # Permission granted to distribute freely for non-commercial purposes only
  7. #
  8. # Copyright (c) 1991, 1992 Jos van der Woude, jvdwoude@hut.nl
  9.  
  10. # If you don't have the gamma() and/or lgamma() functions in your library,
  11. # you can use the following recursive definitions. They are correct for all
  12. # values i / 2 with i = 1, 2, 3, ... This is sufficient for most statistical
  13. # needs.
  14. #logsqrtpi = log(sqrt(pi))
  15. #lgamma(x) = (x<=0.5)?logsqrtpi:((x==1)?0:log(x-1)+lgamma(x-1))
  16. #gamma(x) = exp(lgamma(x))
  17.  
  18. # If you have the lgamma() function compiled into gnuplot, you can use
  19. # alternate definitions for some PDFs. For larger arguments this will result
  20. # in more efficient evalution. Just uncomment the definitions containing the
  21. # string `lgamma', while at the same time commenting out the originals.
  22. # NOTE: In these cases the recursive definition for lgamma() is NOT sufficient!
  23.  
  24. # Some PDFs have alternate definitions of a recursive nature. I suppose these
  25. # are not really very efficient, but I find them aesthetically pleasing to the
  26. # brain.
  27.  
  28. # Define useful constants
  29. fourinvsqrtpi=4.0/sqrt(pi)
  30. invpi=1.0/pi
  31. invsqrt2pi=1.0/sqrt(2.0*pi)
  32. log2=log(2.0)
  33. sqrt2=sqrt(2.0)
  34. sqrt2invpi=sqrt(2.0/pi)
  35. twopi=2.0*pi
  36.  
  37. # define variables plus default values used as parameters in PDFs
  38. # some are integers, others MUST be reals
  39. a=1.0
  40. alpha=0.5
  41. b=2.0
  42. df1=1
  43. df2=1
  44. g=1.0
  45. lambda=1.0
  46. m=0.0
  47. mm=1
  48. mu=0.0
  49. nn=2
  50. n=2
  51. p=0.5
  52. q=0.5
  53. r=1
  54. rho=1.0
  55. sigma=1.0
  56. u=1.0
  57.  
  58. #
  59. #define 1.0/Beta function
  60. #
  61. Binv(p,q)=exp(lgamma(p+q)-lgamma(p)-lgamma(q))
  62.  
  63. #
  64. #define Probability Density Functions (PDFs)
  65. #
  66.  
  67. # NOTE:
  68. # The discrete PDFs are calulated for all real values, using the int()
  69. # function to truncate to integers. This is a monumental waste of processing
  70. # power, but I see no other easy solution. If anyone has any smart ideas
  71. # about this, I would like to know. Setting the sample size to a larger value
  72. # makes the discrete PDFs look better, but takes even more time.
  73.  
  74. # Arcsin PDF
  75. arcsin(x)=invpi/sqrt(r*r-x*x)
  76.  
  77. # Beta PDF
  78. beta(x)=Binv(p,q)*x**(p-1.0)*(1.0-x)**(q-1.0)
  79.  
  80. # Binomial PDF
  81. #binom(x)=n!/(n-int(x))!/int(x)!*p**int(x)*(1.0-p)**(n-int(x))
  82.  
  83. bin_s(x)=n!/(n-int(x))!/int(x)!*p**int(x)*(1.0-p)**(n-int(x))
  84. bin_l(x)=exp(lgamma(n+1)-lgamma(n-int(x)+1)-lgamma(int(x)+1)\
  85. +int(x)*log(p)+(n-int(x))*log(1.0-p))
  86. binom(x)=(n<50)?bin_s(x):bin_l(x)
  87.  
  88. # Cauchy PDF
  89. cauchy(x)=b/(pi*(b*b+(x-a)**2))
  90.  
  91. # Chi-square PDF
  92. #chi(x)=x**(0.5*df1-1.0)*exp(-0.5*x)/gamma(0.5*df1)/2**(0.5*df1)
  93. chi(x)=exp((0.5*df1-1.0)*log(x)-0.5*x-lgamma(0.5*df1)-df1*0.5*log2)
  94.  
  95. # Erlang PDF
  96. erlang(x)=lambda**n/(n-1)!*x**(n-1)*exp(-lambda*x)
  97.  
  98. # Extreme (Gumbel extreme value) PDF
  99. extreme(x)=alpha*(exp(-alpha*(x-u)-exp(-alpha*(x-u))))
  100.  
  101. # F PDF
  102. f(x)=Binv(0.5*df1,0.5*df2)*(df1/df2)**(0.5*df1)*x**(0.5*df1-1.0)/\
  103. (1.0+df1/df2*x)**(0.5*(df1+df2))
  104.  
  105. # Gamma PDF
  106. #g(x)=lambda**rho*x**(rho-1.0)*exp(-lambda*x)/gamma(rho)
  107. g(x)=exp(rho*log(lambda)+(rho-1.0)*log(x)-lgamma(rho)-lambda*x)
  108.  
  109. # Geometric PDF
  110. #geometric(x)=p*(1.0-p)**int(x)
  111. geometric(x)=exp(log(p)+int(x)*log(1.0-p))
  112.  
  113. # Half normal PDF
  114. halfnormal(x)=sqrt2invpi/sigma*exp(-0.5*(x/sigma)**2)
  115.  
  116. # Hypergeometric PDF
  117. hypgeo(x)=(int(x)>mm||int(x)<mm+n-nn)?0:\
  118. mm!/(mm-int(x))!/int(x)!*(nn-mm)!/(n-int(x))!/(nn-mm-n+int(x))!*(nn-n)!*n!/nn!
  119.  
  120. # Laplace PDF
  121. laplace(x)=0.5/b*exp(-abs(x-a)/b)
  122.  
  123. # Logistic PDF
  124. logistic(x)=lambda*exp(-lambda*(x-a))/(1.0+exp(-lambda*(x-a)))**2
  125.  
  126. # Lognormal PDF
  127. lognormal(x)=invsqrt2pi/sigma/x*exp(-0.5*((log(x)-mu)/sigma)**2)
  128.  
  129. # Maxwell PDF
  130. maxwell(x)=fourinvsqrtpi*a**3*x*x*exp(-a*a*x*x)
  131.  
  132. # Negative binomial PDF
  133. #negbin(x)=(r+int(x)-1)!/int(x)!/(r-1)!*p**r*(1.0-p)**int(x)
  134. negbin(x)=exp(lgamma(r+int(x))-lgamma(r)-lgamma(int(x)+1)+\
  135. r*log(p)+int(x)*log(1.0-p))
  136.  
  137. # Negative exponential PDF
  138. nexp(x)=lambda*exp(-lambda*x)
  139.  
  140. # Normal PDF
  141. normal(x)=invsqrt2pi/sigma*exp(-0.5*((x-mu)/sigma)**2)
  142.  
  143. # Pareto PDF
  144. pareto(x)=x<a?0:b/x*(a/x)**b
  145.  
  146. # Poisson PDF
  147. poisson(x)=mu**int(x)/int(x)!*exp(-mu)
  148. #poisson(x)=exp(int(x)*log(mu)-lgamma(int(x)+1)-mu)
  149. #poisson(x)=(x<1)?exp(-mu):mu/int(x)*poisson(x-1)
  150. #lpoisson(x)=(x<1)?-mu:log(mu)-log(int(x))+lpoisson(x-1)
  151.  
  152. # Rayleigh PDF
  153. rayleigh(x)=lambda*2.0*x*exp(-lambda*x*x)
  154.  
  155. # Sine PDF
  156. sine(x)=2.0/a*sin(n*pi*x/a)**2
  157.  
  158. # t (Student's t) PDF
  159. t(x)=Binv(0.5*df1,0.5)/sqrt(df1)*(1.0+(x*x)/df1)**(-0.5*(df1+1.0))
  160.  
  161. # Triangular PDF
  162. triangular(x)=1.0/g-abs(x-m)/(g*g)
  163.  
  164. # Uniform PDF
  165. uniform(x)=1.0/(b-a)
  166.  
  167. # Weibull PDF
  168. weibull(x)=lambda*n*x**(n-1)*exp(-lambda*x**n)
  169.  
  170. #
  171. #define Cumulative Distribution Functions (CDFs)
  172. #
  173.  
  174. # Arcsin CDF
  175. carcsin(x)=0.5+invpi*asin(x/r)
  176.  
  177. # incomplete Beta CDF
  178. cbeta(x)=ibeta(p,q,x)
  179.  
  180. # Binomial CDF
  181. #cbinom(x)=(x<1)?binom(0):binom(x)+cbinom(x-1)
  182. cbinom(x)=ibeta(n-x,x+1.0,1.0-p)
  183.  
  184. # Cauchy CDF
  185. ccauchy(x)=0.5+invpi*atan((x-a)/b)
  186.  
  187. # Chi-square CDF
  188. cchi(x)=igamma(0.5*df1,0.5*x)
  189.  
  190. # Erlang CDF
  191. # approximation, using first three terms of expansion
  192. cerlang(x)=1.0-exp(-lambda*x)*(1.0+lambda*x+0.5*(lambda*x)**2)
  193.  
  194. # Extreme (Gumbel extreme value) CDF
  195. cextreme(x)=exp(-exp(-alpha*(x-u)))
  196.  
  197. # F CDF
  198. cf(x)=1.0-ibeta(0.5*df2,0.5*df1,df2/(df2+df1*x))
  199.  
  200. # incomplete Gamma CDF
  201. cgamma(x)=igamma(rho,x)
  202.  
  203. # Geometric CDF
  204. cgeometric(x)=(x<1)?geometric(0):geometric(x)+cgeometric(x-1)
  205.  
  206. # Half normal CDF
  207. chalfnormal(x)=erf(x/sigma/sqrt2)
  208.  
  209. # Hypergeometric CDF
  210. chypgeo(x)=(x<1)?hypgeo(0):hypgeo(x)+chypgeo(x-1)
  211.  
  212. # Laplace CDF
  213. claplace(x)=(x<a)?0.5*exp((x-a)/b):1.0-0.5*exp(-(x-a)/b)
  214.  
  215. # Logistic CDF
  216. clogistic(x)=1.0/(1.0+exp(-lambda*(x-a)))
  217.  
  218. # Lognormal CDF
  219. clognormal(x)=cnormal(log(x))
  220.  
  221. # Maxwell CDF
  222. cmaxwell(x)=igamma(1.5,a*a*x*x)
  223.  
  224. # Negative binomial CDF
  225. cnegbin(x)=(x<1)?negbin(0):negbin(x)+cnegbin(x-1)
  226.  
  227. # Negative exponential CDF
  228. cnexp(x)=1.0-exp(-lambda*x)
  229.  
  230. # Normal CDF
  231. cnormal(x)=0.5+0.5*erf((x-mu)/sigma/sqrt2)
  232. #cnormal(x)=0.5+((x>mu)?0.5:-0.5)*igamma(0.5,0.5*((x-mu)/sigma)**2)
  233.  
  234. # Pareto CDF
  235. cpareto(x)=x<a?0:1.0-(a/x)**b
  236.  
  237. # Poisson CDF
  238. #cpoisson(x)=(x<1)?poisson(0):poisson(x)+cpoisson(x-1)
  239. cpoisson(x)=1.0-igamma(x+1.0,mu)
  240.  
  241. # Rayleigh CDF
  242. crayleigh(x)=1.0-exp(-lambda*x*x)
  243.  
  244. # Sine CDF
  245. csine(x)=x/a-sin(n*twopi*x/a)/(n*twopi)
  246.  
  247. # t (Student's t) CDF
  248. ct(x)=(x<0.0)?0.5*ibeta(0.5*df1,0.5,df1/(df1+x*x)):\
  249. 1.0-0.5*ibeta(0.5*df1,0.5,df1/(df1+x*x))
  250.  
  251. # Triangular PDF
  252. ctriangular(x)=0.5+(x-m)/g-(x-m)*abs(x-m)/(2.0*g*g)
  253.  
  254. # Uniform CDF
  255. cuniform(x)=(x-a)/(b-a)
  256.  
  257. # Weibull CDF
  258. cweibull(x)=1.0-exp(-lambda*x**n)
  259.